####################
# Loading rep2 data
####################

ldat_new <- read.loom.matrices("/Volumes/bioinfomatics$/Mehdi/Matthias/scRNA-seq/scRNAseq_CRG_02012020/onefilepercell_HFWL7BBXY_1_IDT-DUI-NXT-66_and_others_0VWWZ.loom")
## reading loom file via hdf5r...
ldat_new <- lapply(ldat_new,function(x) {
  colnames(x) <-  gsub(".bam","",gsub("onefilepercell_HFWL7BBXY_1_IDT-DUI-NXT-66_and_others_0VWWZ:","",colnames(x)))
  x
})

meta_data <- read.csv(file="/Volumes/bioinfomatics$/Mehdi/Matthias/scRNA-seq/scRNAseq_CRG_02012020/sample_info.csv", header=TRUE, sep=",")
rownames(meta_data) <- meta_data$sample_ID
meta_data$rep <- 'rep2'
meta_data <- meta_data[,c('rep','sample_label')]
colnames(meta_data) <- c('rep','celltype')
meta_data$celltype <- gsub('TCRint_DP_WT','TCRhiDP',gsub('CD8_SP_TCRpos_WT','CD8SP',gsub("CD4_SP_WT","CD4SP",gsub("CD4_pos_8lo_WT","CD4+8low",gsub("CD69_pos_DP_WT","CD69posDP",gsub("CD69_neg_DP_WT","CD69negDP",meta_data$celltype))))))


meta_data$celltype <- gsub('TCRint_DP_MHCclassII_KO','TCRhiDP_C2KO',gsub('CD8_SP_TCRpos_MHCclassII_KO','CD8SP_C2KO',gsub("CD4_SP_MHCclassII_KO","CD4SP_C2KO",gsub("CD4_pos_8lo_classII_KO","CD4+8low_C2KO",gsub("CD69_pos_DP_MHCclassII_KO","CD69posDP_C2KO",gsub("CD69_neg_DP_MHCclassII_KO","CD69negDP_C2KO",meta_data$celltype))))))

meta_data <- meta_data[!(rownames(meta_data)=="HFWL7BBXY_1_IDT-DUI-NXT-49" | rownames(meta_data)=="HFWL7BBXY_1_IDT-DUI-NXT-333"),]
meta_data[grep('C2KO',meta_data$celltype),]$rep <- 'rep2_KO'
meta_data$celltype <- gsub('_C2KO','', meta_data$celltype)

# exonic read (spliced) expression matrix
emat_rep2 <- ldat_new$spliced
emat_rep2 <- emat_rep2[, colnames(emat_rep2) %in% as.vector(rownames(meta_data))];


# intronic read (unspliced) expression matrix
nmat_rep2 <- ldat_new$unspliced
nmat_rep2 <- nmat_rep2[, colnames(nmat_rep2) %in% as.vector(rownames(meta_data))];

# spanning read (intron+exon) expression matrix
smat_rep2 <- ldat_new$spanning;
smat_rep2 <- smat_rep2[, colnames(smat_rep2) %in% as.vector(rownames(meta_data))];


######################
# Loading rep1 data
######################

ldat <- read.loom.matrices("/Volumes/bioinfomatics$/Mehdi/Matthias/scRNA-seq/onefilepercell_P2221_N707-S505_and_others_96WOS.loom")
## reading loom file via hdf5r...
ldat <- lapply(ldat,function(x) {
  colnames(x) <-  gsub(".bam","",gsub(".*:","",colnames(x)))
  x
})

WishBone_tsne <- read.csv(file="/Volumes/bioinfomatics$/Mehdi/Matthias/scRNA-seq/WishBone_tsne_with_highTCR.csv", header=TRUE, sep=",")
colnames(WishBone_tsne) <- c('CellName','x','y','Sample')

colData <- read.csv(file="/Volumes/bioinfomatics$/Mehdi/Matthias/scRNA-seq/colData.csv", header=TRUE, sep=",")

merged_meta1 <- merge(x = WishBone_tsne, y = colData, by = "CellName")
color_dt <- data.frame(Sample.x=c('WT1','WT2','WT3','WT4','WT5','highTCR_Rest'), color=c('#f44e42','#ffac38','#bcff37','#36eeff','#9335ff','#36ffee'))
merged_meta <- merge(x = merged_meta1, y = color_dt, by = "Sample.x")
merged_meta <- merged_meta[,c('CellCode','Sample.x')]
rownames(merged_meta) <- merged_meta$CellCode
merged_meta$rep <- 'rep1'
merged_meta <- merged_meta[,c('rep','Sample.x')]
colnames(merged_meta) <- c('rep','celltype')
merged_meta$celltype <- gsub('highTCR_Rest','TCRhiDP',gsub('WT5','CD8SP',gsub("WT4","CD4SP",gsub("WT3","CD4+8low",gsub("WT2","CD69posDP",gsub("WT1","CD69negDP",merged_meta$celltype))))))

# exonic read (spliced) expression matrix
emat_rep1 <- ldat$spliced
colnames(emat_rep1) <- gsub("-","_",colnames(emat_rep1))
emat_rep1 <- emat_rep1[, colnames(emat_rep1) %in% as.vector(rownames(merged_meta))];


# intronic read (unspliced) expression matrix
nmat_rep1 <- ldat$unspliced
colnames(nmat_rep1) <- gsub("-","_",colnames(nmat_rep1))
nmat_rep1 <- nmat_rep1[, colnames(nmat_rep1) %in% as.vector(rownames(merged_meta))];

# spanning read (intron+exon) expression matrix
smat_rep1 <- ldat$spanning;
colnames(smat_rep1) <- gsub("-","_",colnames(smat_rep1))
smat_rep1 <- smat_rep1[, colnames(smat_rep1) %in% as.vector(rownames(merged_meta))];


############
#merging rep1 & rep2
############

emat <- cbind(emat_rep1,emat_rep2)
nmat <- cbind(nmat_rep1,nmat_rep2)
smat <- cbind(smat_rep1,smat_rep2)
meta_info <- rbind(merged_meta,meta_data)


###########
# Running Seurat 
###########

WT_Thymocye <- CreateSeuratObject(emat, meta.data = meta_info)
WT_Thymocye.list <- SplitObject(WT_Thymocye, split.by = "rep")

for (i in 1:length(WT_Thymocye.list)) {
    WT_Thymocye.list[[i]] <- NormalizeData(WT_Thymocye.list[[i]], verbose = FALSE)
    WT_Thymocye.list[[i]] <- FindVariableFeatures(WT_Thymocye.list[[i]], selection.method = "vst", 
        nfeatures = 2000, verbose = FALSE)
}

Rep1_log <- t(GetAssayData(WT_Thymocye.list[[1]],slot='data',assay = "RNA"))
Rep2_log <- t(GetAssayData(WT_Thymocye.list[[3]],slot='data',assay = "RNA"))
Rep2KO_log <- t(GetAssayData(WT_Thymocye.list[[2]],slot='data',assay = "RNA"))

meta_info$x1 <- 'Others'
meta_info$x2 <- 'Others'
meta_info$x3 <- 'Others'
meta_info$x4 <- 'Others'
meta_info$x5 <- 'Cd8-'
meta_info$x6 <- 'Cd4-'
meta_info$x7 <- 'Itm2a-'
meta_info$x8 <- 'Stat1-'


colnames(meta_info) <- c("rep", "celltype", "Cd4+Cd8+","Cd4+Cd8-", "Cd4-Cd8+", "Cd4-Cd8-", 'Cd8', 'Cd4', 'Itm2a', 'Stat1') 

meta_info[rownames(Rep1_log[Rep1_log[,'Cd4'] > 0.1 & Rep1_log[,'Cd8a'] > 0.1,]),'Cd4+Cd8+'] <- 'Cd4+Cd8+'
meta_info[rownames(Rep2_log[Rep2_log[,'Cd4'] > 0.1 & Rep2_log[,'Cd8a'] > 0.1,]),'Cd4+Cd8+'] <- 'Cd4+Cd8+'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Cd4'] > 1 & Rep2KO_log[,'Cd8a'] > 1,]),'Cd4+Cd8+'] <- 'Cd4+Cd8+'

meta_info[rownames(Rep1_log[Rep1_log[,'Cd4'] > 0.1 & Rep1_log[,'Cd8a'] < 0.1,]),'Cd4+Cd8-'] <- 'Cd4+Cd8-'
meta_info[rownames(Rep2_log[Rep2_log[,'Cd4'] > 0.1 & Rep2_log[,'Cd8a'] < 0.1,]),'Cd4+Cd8-'] <- 'Cd4+Cd8-'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Cd4'] > 1 & Rep2KO_log[,'Cd8a'] < 1,]),'Cd4+Cd8-'] <- 'Cd4+Cd8-'

meta_info[rownames(Rep1_log[Rep1_log[,'Cd4'] < 0.1 & Rep1_log[,'Cd8a'] > 0.1,]),'Cd4-Cd8+'] <- 'Cd4-Cd8+'
meta_info[rownames(Rep2_log[Rep2_log[,'Cd4'] < 0.1 & Rep2_log[,'Cd8a'] > 0.1,]),'Cd4-Cd8+'] <- 'Cd4-Cd8+'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Cd4'] < 1 & Rep2KO_log[,'Cd8a'] > 1,]),'Cd4-Cd8+'] <- 'Cd4-Cd8+'

meta_info[rownames(Rep1_log[Rep1_log[,'Cd4'] < 0.1 & Rep1_log[,'Cd8a'] < 0.1,]),'Cd4-Cd8-'] <- 'Cd4-Cd8-'
meta_info[rownames(Rep2_log[Rep2_log[,'Cd4'] < 0.1 & Rep2_log[,'Cd8a'] < 0.1,]),'Cd4-Cd8-'] <- 'Cd4-Cd8-'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Cd4'] < 1 & Rep2KO_log[,'Cd8a'] < 1,]),'Cd4-Cd8-'] <- 'Cd4-Cd8-'

meta_info[rownames(Rep1_log[Rep1_log[,'Cd8a'] > 0.1,]),'Cd8'] <- 'Cd8+'
meta_info[rownames(Rep2_log[Rep2_log[,'Cd8a'] > 0.1,]),'Cd8'] <- 'Cd8+'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Cd8a'] > 0.1,]),'Cd8'] <- 'Cd8+'

meta_info[rownames(Rep1_log[Rep1_log[,'Cd4'] > 0.1,]),'Cd4'] <- 'Cd4+'
meta_info[rownames(Rep2_log[Rep2_log[,'Cd4'] > 0.1,]),'Cd4'] <- 'Cd4+'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Cd4'] > 0.1,]),'Cd4'] <- 'Cd4+'

meta_info[rownames(Rep1_log[Rep1_log[,'Itm2a'] > 0.5,]),'Itm2a'] <- 'Itm2a+'
meta_info[rownames(Rep2_log[Rep2_log[,'Itm2a'] > 0.5,]),'Itm2a'] <- 'Itm2a+'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Itm2a'] > 0.1,]),'Itm2a'] <- 'Itm2a+'

meta_info[rownames(Rep1_log[Rep1_log[,'Stat1'] > 0.2,]),'Stat1'] <- 'Stat1+'
meta_info[rownames(Rep2_log[Rep2_log[,'Stat1'] > 0.2,]),'Stat1'] <- 'Stat1+'
meta_info[rownames(Rep2KO_log[Rep2KO_log[,'Stat1'] > 0.1,]),'Stat1'] <- 'Stat1+'

meta_info[meta_info$celltype=='CD4SP' & meta_info$Itm2a=='Itm2a+' & meta_info$Stat1=='Stat1-', 'celltype'] <- 'CD4SP_Immature'
meta_info[meta_info$celltype=='CD4SP' & (meta_info$Itm2a=='Itm2a-' | meta_info$Stat1=='Stat1+'), 'celltype'] <- 'CD4SP_Mature'



WT_Thymocye <- CreateSeuratObject(emat, meta.data = meta_info)
WT_Thymocye.list <- SplitObject(WT_Thymocye, split.by = "rep")

for (i in 1:length(WT_Thymocye.list)) {
    WT_Thymocye.list[[i]] <- NormalizeData(WT_Thymocye.list[[i]], verbose = FALSE)
    WT_Thymocye.list[[i]] <- FindVariableFeatures(WT_Thymocye.list[[i]], selection.method = "vst", 
        nfeatures = 2000, verbose = FALSE)
}

meta_info <- WT_Thymocye@meta.data

reference.list <- WT_Thymocye.list[c("rep1", "rep2","rep2_KO")]
WT_Thymocye.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

WT_Thymocye.integrated <- IntegrateData(anchorset = WT_Thymocye.anchors, dims = 1:30)

DefaultAssay(WT_Thymocye.integrated) <- "integrated"

# Run the standard workflow for visualization and clustering
WT_Thymocye.integrated <- ScaleData(WT_Thymocye.integrated, verbose = FALSE)
WT_Thymocye.integrated <- RunPCA(WT_Thymocye.integrated, npcs = 30, verbose = FALSE)
WT_Thymocye.integrated <- RunUMAP(WT_Thymocye.integrated, reduction = "pca", dims = 1:30)

gene_set <- read.csv(file="/Volumes/bioinfomatics$/Mehdi/Matthias/scRNA-seq/scRNAseq_CRG_02012020/gene_set/DE_gene_drodeRes_rep2_fset.csv", header=FALSE, sep=",")

emat <- cbind(emat_rep1,emat_rep2)
nmat <- cbind(nmat_rep1,nmat_rep2)
smat <- cbind(smat_rep1,smat_rep2)
# meta_info <- rbind(merged_meta,meta_data)

emat <- emat[intersect(gene_set$V1,rownames(emat)),]
emat <- emat[!(rownames(emat) =='Cd4' | rownames(emat) == 'Cd8a' | rownames(emat) =='Cd8b1'),]

###########
#Running Seurat 
###########

WT_Thymocye <- CreateSeuratObject(emat, meta.data = meta_info)
WT_Thymocye.list <- SplitObject(WT_Thymocye, split.by = "rep")

for (i in 1:length(WT_Thymocye.list)) {
    WT_Thymocye.list[[i]] <- NormalizeData(WT_Thymocye.list[[i]], verbose = FALSE)
    WT_Thymocye.list[[i]] <- FindVariableFeatures(WT_Thymocye.list[[i]], selection.method = "vst", 
        nfeatures = 2000, verbose = FALSE)
}

reference.list <- WT_Thymocye.list[c("rep1", "rep2","rep2_KO")]
WT_Thymocye.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

WT_Thymocye.integrated <- IntegrateData(anchorset = WT_Thymocye.anchors, dims = 1:30)

# switch to integrated assay. The variable features of this assay are automatically
# set during IntegrateData

DefaultAssay(WT_Thymocye.integrated) <- "integrated"

# Run the standard workflow for visualization and clustering
WT_Thymocye.integrated <- ScaleData(WT_Thymocye.integrated, verbose = FALSE)
WT_Thymocye.integrated <- RunPCA(WT_Thymocye.integrated, npcs = 30, verbose = FALSE)
WT_Thymocye.integrated <- RunUMAP(WT_Thymocye.integrated, reduction = "pca", dims = 1:30)

# Filtering cells with too many or too few cells

Rep1_filtered_rep2 <- subset(WT_Thymocye.integrated, cells = rownames(meta_info[meta_info$rep=='rep1' |
         (meta_info$rep=='rep2_KO' & meta_info$celltype=='CD69posDP' & meta_info$nCount_RNA < 1500000 & meta_info$nCount_RNA > 750000) |
         (meta_info$rep=='rep2_KO' & meta_info$celltype!='CD69posDP' & meta_info$nCount_RNA < 1500000 & meta_info$nCount_RNA > 500000) |
         (meta_info$rep=='rep2' & meta_info$celltype=='CD69posDP' & meta_info$nCount_RNA < 1500000 & meta_info$nCount_RNA > 750000) |
         (meta_info$rep=='rep2' & meta_info$celltype!='CD69posDP' & meta_info$nCount_RNA < 1500000 & meta_info$nCount_RNA > 500000) 
                ,]))


# Dividing CD4SP to CD4SP Immature and Mature

meta_info_new <- Rep1_filtered_rep2@meta.data
PCA <- as.data.frame(Rep1_filtered_rep2@reductions$pca@cell.embeddings)
rownames(PCA) <-  rownames(Rep1_filtered_rep2@reductions$pca@cell.embeddings)

meta_info_new$PC_1 <- PCA$PC_1
meta_info_new$PC_2 <- PCA$PC_2
meta_info_new$PC_3 <- PCA$PC_3

meta_info_new[meta_info_new$celltype=='CD4+8low' & meta_info_new$PC_2 < meta_info_new$PC_1,'celltype'] <- 'CD4SP_Immature'

meta_info[rownames(meta_info_new),'celltype'] <- meta_info_new$celltype
gene_set <- read.csv(file="/Volumes/bioinfomatics$/Mehdi/Matthias/scRNA-seq/scRNAseq_CRG_02012020/gene_set/DE_gene_drodeRes_rep2_fset.csv", header=FALSE, sep=",")

emat <- cbind(emat_rep1,emat_rep2)
nmat <- cbind(nmat_rep1,nmat_rep2)
smat <- cbind(smat_rep1,smat_rep2)
# meta_info <- rbind(merged_meta,meta_data)

emat <- emat[intersect(gene_set$V1,rownames(emat)),]
emat <- emat[!(rownames(emat) =='Cd4' | rownames(emat) == 'Cd8a' | rownames(emat) =='Cd8b1'),]

# gene_set <- gene_set[!(gene_set$V1 =='Cd4' | gene_set$V1 =='Cd8a' | gene_set$V1 =='Cd8b1'),] 

###########
#Running Seurat 
###########

WT_Thymocye <- CreateSeuratObject(emat, meta.data = meta_info)
WT_Thymocye.list <- SplitObject(WT_Thymocye, split.by = "rep")

for (i in 1:length(WT_Thymocye.list)) {
    WT_Thymocye.list[[i]] <- NormalizeData(WT_Thymocye.list[[i]], verbose = FALSE)
    WT_Thymocye.list[[i]] <- FindVariableFeatures(WT_Thymocye.list[[i]], selection.method = "vst", 
        nfeatures = 2000, verbose = FALSE)
}

reference.list <- WT_Thymocye.list[c("rep1", "rep2", "rep2_KO")]
WT_Thymocye.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

WT_Thymocye.integrated <- IntegrateData(anchorset = WT_Thymocye.anchors, dims = 1:30)


library(ggplot2)
library(cowplot)
# switch to integrated assay. The variable features of this assay are automatically
# set during IntegrateData
DefaultAssay(WT_Thymocye.integrated) <- "integrated"

# Run the standard workflow for visualization and clustering
WT_Thymocye.integrated <- ScaleData(WT_Thymocye.integrated, verbose = FALSE)
WT_Thymocye.integrated <- RunPCA(WT_Thymocye.integrated, npcs = 30, verbose = FALSE)
WT_Thymocye.integrated <- RunUMAP(WT_Thymocye.integrated, reduction = "pca", dims = 1:30)

PCA <- as.data.frame(WT_Thymocye.integrated@reductions$pca@cell.embeddings)



# Original runing 
emat <- cbind(emat_rep1,emat_rep2)
nmat <- cbind(nmat_rep1,nmat_rep2)
smat <- cbind(smat_rep1,smat_rep2)

###########
#Running Seurat 
###########

WT_Thymocye <- CreateSeuratObject(emat, meta.data = meta_info)
WT_Thymocye.list <- SplitObject(WT_Thymocye, split.by = "rep")

for (i in 1:length(WT_Thymocye.list)) {
    WT_Thymocye.list[[i]] <- NormalizeData(WT_Thymocye.list[[i]], verbose = FALSE)
    WT_Thymocye.list[[i]] <- FindVariableFeatures(WT_Thymocye.list[[i]], selection.method = "vst", 
        nfeatures = 5000, verbose = FALSE)
}

reference.list <- WT_Thymocye.list[c("rep1", "rep2", "rep2_KO")]
WT_Thymocye.anchors <- FindIntegrationAnchors(object.list = reference.list, dims = 1:30)

WT_Thymocye.integrated <- IntegrateData(anchorset = WT_Thymocye.anchors, dims = 1:30)


library(ggplot2)
library(cowplot)
# switch to integrated assay. The variable features of this assay are automatically
# set during IntegrateData
DefaultAssay(WT_Thymocye.integrated) <- "integrated"

WT_Thymocye.integrated <- ScaleData(WT_Thymocye.integrated, verbose = FALSE)
WT_Thymocye.integrated <- RunPCA(WT_Thymocye.integrated, npcs = 30, verbose = FALSE)
WT_Thymocye.integrated <- RunUMAP(WT_Thymocye.integrated, reduction = "pca", dims = 1:30)


WT_Thymocye.integrated@reductions$pca@cell.embeddings <-  as.matrix(PCA)


Rep1_filtered_rep2 <- subset(WT_Thymocye.integrated, cells = rownames(meta_info[
         (meta_info$rep=='rep2' & meta_info$celltype=='CD69posDP' & meta_info$nCount_RNA < 1500000 & meta_info$nCount_RNA > 750000) |
         (meta_info$rep=='rep2' & meta_info$celltype!='CD69posDP' & meta_info$nCount_RNA < 1500000 & meta_info$nCount_RNA > 500000) 
                ,]))

DefaultAssay(Rep1_filtered_rep2) <- "RNA"
Rep1_filtered_rep2@meta.data$CD4_expression <- 0
Rep1_filtered_rep2@meta.data$CD4_expression <- t(GetAssayData(Rep1_filtered_rep2,slot='data',assay = "RNA"))[,"Cd4"]

Rep1_filtered_rep2@meta.data$Itm2a_expression <- 0
Rep1_filtered_rep2@meta.data$Itm2a_expression <- t(GetAssayData(Rep1_filtered_rep2,slot='data',assay = "RNA"))[,"Itm2a"]

Rep1_filtered_rep2@meta.data$Tbx21_expression <- 0
Rep1_filtered_rep2@meta.data$Tbx21_expression <- t(GetAssayData(Rep1_filtered_rep2,slot='data',assay = "RNA"))[,"Tbx21"]

draw_FeaturePlot <- function(pair) {
P1 <- FeaturePlot(Rep1_filtered_rep2,features = pair, slot='data',label = F, blend = TRUE, order = TRUE, reduction = "pca", pt.size=1, cols = c("darkgrey", "#ff0000", "#00ff00"), blend.threshold = 0.1)
P2 <- FeaturePlot(Rep1_filtered_rep2,features = pair, slot='data',label = F, blend = TRUE, order = TRUE, reduction = "pca", dims = c(1,3), pt.size=1, cols = c("darkgrey", "#ff0000", "#00ff00"), blend.threshold = 0.1)
print(P1)
print(P2)
}

0.1 Figure 2: scRNA-seq captures maturation, activation, and lineage identity as principal components of CD4/CD8 lineage choice and differentiation

draw_FeaturePlot(c("Rag1","Ccr7"))

draw_FeaturePlot(c("Cd69","Cd69"))

draw_FeaturePlot(c("CD4_expression","Cd8a"))

0.2 Figure 6.ab) CD4 and CD8 lineage characteristics are established at different times during CD4/CD8 lineage choice and differentiation

draw_FeaturePlot(c("Itm2a_expression","Zbtb7b"))

draw_FeaturePlot(c("Itm2a_expression","Cd40lg"))

draw_FeaturePlot(c("Itm2a_expression","Nkg7"))

draw_FeaturePlot(c("Itm2a_expression","Itgae"))

0.3 Sup Figure 2.bcde) Key genes and correlations

draw_FeaturePlot(c("Rag1","Dntt"))

draw_FeaturePlot(c("Rag1","Cd69"))

draw_FeaturePlot(c("Ccr7","Il7r"))

draw_FeaturePlot(c("CD4_expression","Cd8a"))

0.4 Session Info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
##  [1] grid      splines   parallel  stats4    stats     graphics  grDevices
##  [8] utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.1.0            gridExtra_2.3              
##  [3] cowplot_1.0.0               tiff_0.1-5                 
##  [5] jpeg_0.1-8.1                VennDiagram_1.6.20         
##  [7] futile.logger_1.4.3         metaMA_3.1.2               
##  [9] SCORPIUS_1.0.6              dynamicTreeCut_1.63-1      
## [11] cluster_2.1.0               org.Mm.eg.db_3.10.0        
## [13] AnnotationDbi_1.48.0        knitr_1.28                 
## [15] scran_1.14.6                scater_1.14.6              
## [17] RColorBrewer_1.1-2          Seurat_3.1.4               
## [19] DT_0.12                     slingshot_1.4.0            
## [21] princurve_2.1.4             splatter_1.10.1            
## [23] SingleCellExperiment_1.8.0  gam_1.16.1                 
## [25] foreach_1.4.8               ggpubr_0.2.5               
## [27] magrittr_1.5                DESeq2_1.26.0              
## [29] SummarizedExperiment_1.16.1 DelayedArray_0.12.2        
## [31] BiocParallel_1.20.1         matrixStats_0.56.0         
## [33] Biobase_2.46.0              GenomicRanges_1.38.0       
## [35] GenomeInfoDb_1.22.0         IRanges_2.20.2             
## [37] S4Vectors_0.24.3            BiocGenerics_0.32.0        
## [39] velocyto.R_0.6              Matrix_1.2-18              
## [41] forcats_0.5.0               stringr_1.4.0              
## [43] dplyr_0.8.5                 purrr_0.3.3                
## [45] readr_1.3.1                 tidyr_1.0.2                
## [47] tibble_2.1.3                ggplot2_3.3.0              
## [49] tidyverse_1.3.0             dyno_0.1.2                 
## [51] dynwrap_1.2.0.9000          dynplot_1.0.2.9000         
## [53] dynmethods_1.0.5            dynguidelines_1.0.0        
## [55] dynfeature_1.0.0.9000      
## 
## loaded via a namespace (and not attached):
##   [1] rsvd_1.0.3               Hmisc_4.3-1              ica_1.0-2               
##   [4] ps_1.3.2                 lmtest_0.9-37            crayon_1.3.4            
##   [7] MASS_7.3-51.5            nlme_3.1-144             backports_1.1.5         
##  [10] reprex_0.3.0             rlang_0.4.5              GA_3.2                  
##  [13] XVector_0.26.0           ROCR_1.0-7               readxl_1.3.1            
##  [16] irlba_2.3.3              limma_3.42.2             dyndimred_1.0.3.9000    
##  [19] bit64_0.9-7              glue_1.3.2               pheatmap_1.0.12         
##  [22] sctransform_0.2.1        processx_3.4.2           vipor_0.4.5             
##  [25] dynutils_1.0.5           haven_2.2.0              tidyselect_1.0.0        
##  [28] lmds_0.1.0               fitdistrplus_1.0-14      XML_3.99-0.3            
##  [31] zoo_1.8-7                xtable_1.8-4             rje_1.10.15             
##  [34] evaluate_0.14            bibtex_0.4.2.2           Rdpack_0.11-1           
##  [37] cli_2.0.2                zlibbioc_1.32.0          sn_1.5-5                
##  [40] rstudioapi_0.11          rpart_4.1-15             lambda.r_1.2.4          
##  [43] shiny_1.4.0.2            BiocSingular_1.2.2       xfun_0.12               
##  [46] multtest_2.42.0          caTools_1.18.0           tidygraph_1.1.2         
##  [49] TSP_1.1-9                pcaMethods_1.78.0        ggrepel_0.8.2           
##  [52] ape_5.3                  listenv_0.8.0            png_0.1-7               
##  [55] future_1.16.0            withr_2.1.2              bitops_1.0-6            
##  [58] ggforce_0.3.1            ranger_0.12.1            plyr_1.8.6              
##  [61] cellranger_1.1.0         dqrng_0.2.1              pillar_1.4.3            
##  [64] RcppParallel_5.0.0       gplots_3.0.3             multcomp_1.4-12         
##  [67] fs_1.3.2                 DelayedMatrixStats_1.8.0 vctrs_0.2.4             
##  [70] generics_0.0.2           metap_1.3                tools_3.6.3             
##  [73] foreign_0.8-75           beeswarm_0.2.3           munsell_0.5.0           
##  [76] tweenr_1.0.1             fastmap_1.0.1            compiler_3.6.3          
##  [79] httpuv_1.5.2             babelwhale_1.0.1         plotly_4.9.2            
##  [82] GenomeInfoDbData_1.2.2   edgeR_3.28.1             lattice_0.20-38         
##  [85] mutoss_0.1-12            later_1.0.0              jsonlite_1.6.1          
##  [88] scales_1.1.0             pbapply_1.4-2            genefilter_1.68.0       
##  [91] lazyeval_0.2.2           promises_1.1.0           latticeExtra_0.6-29     
##  [94] reticulate_1.13          checkmate_2.0.0          rmarkdown_2.1           
##  [97] sandwich_2.5-1           webshot_0.5.2            statmod_1.4.34          
## [100] Rtsne_0.15               uwot_0.1.5               igraph_1.2.4.2          
## [103] proxyC_0.1.5             survival_3.1-8           numDeriv_2016.8-1.1     
## [106] yaml_2.2.1               plotrix_3.7-7            carrier_0.1.0           
## [109] htmltools_0.4.0          memoise_1.1.0            locfit_1.5-9.1          
## [112] graphlayouts_0.6.0       viridisLite_0.3.0        digest_0.6.25           
## [115] assertthat_0.2.1         mime_0.9                 futile.options_1.0.1    
## [118] npsurv_0.4-0             dynparam_1.0.0           RSQLite_2.2.0           
## [121] future.apply_1.4.0       lsei_1.2-0               remotes_2.1.1           
## [124] data.table_1.12.8        blob_1.2.1               labeling_0.3            
## [127] Formula_1.2-3            RCurl_1.98-1.1           SMVar_1.3.3             
## [130] broom_0.5.5              hms_0.5.3                modelr_0.1.6            
## [133] colorspace_1.4-1         base64enc_0.1-3          mnormt_1.5-6            
## [136] ggbeeswarm_0.6.0         nnet_7.3-12              Rcpp_1.0.3              
## [139] mclust_5.4.5             RANN_2.6.1               mvtnorm_1.1-0           
## [142] fansi_0.4.1              R6_2.4.1                 ggridges_0.5.2          
## [145] lifecycle_0.2.0          acepack_1.4.1            formatR_1.7             
## [148] TFisher_0.2.0            ggsignif_0.6.0           gdata_2.18.0            
## [151] leiden_0.3.3             testthat_2.3.2           RcppAnnoy_0.0.16        
## [154] TH.data_1.0-10           iterators_1.0.12         htmlwidgets_1.5.1       
## [157] polyclip_1.10-0          rvest_0.3.5              mgcv_1.8-31             
## [160] globals_0.12.5           htmlTable_1.13.3         patchwork_1.0.0.9000    
## [163] codetools_0.2-16         lubridate_1.7.4          gtools_3.8.1            
## [166] dbplyr_1.4.2             gtable_0.3.0             tsne_0.1-3              
## [169] DBI_1.1.0                httr_1.4.1               KernSmooth_2.23-16      
## [172] stringi_1.4.6            reshape2_1.4.3           farver_2.0.3            
## [175] annotate_1.64.0          viridis_0.5.1            xml2_1.2.5              
## [178] BiocNeighbors_1.4.2      geneplotter_1.64.0       bit_1.1-15.2            
## [181] ggraph_2.0.1             pkgconfig_2.0.3          gbRd_0.4-11